C   GENERAL PLOTTING PACKAGE FOR PRINTER OR TERMINAL
C
      SUBROUTINE PPINIT(NTYPE,IFLU,XMN,XMX,YMN,YMX)
C
C   INITIALIZATION ROUTINE FOR PRINTER PLOTTING PACKAGE
C
C
C      WHERE:
C
C     NTYPE - PLOT TYPE: 1 = LINEAR
C                        2 = SCALED SEMI-LOG (Y LINEAR)
C                        3 = SCALED SEMI-LOG (X LINEAR)
C                        4 = SCALED LOG-LOG
C                        5 = SEMI-LOG (Y LINEAR)
C                        6 = SEMI-LOG (X LINEAR)
C                        7 = LOG-LOG
C      IFLU - FORTRAN LOGICAL UNIT USED FOR PLOT OUTPUT
C       XMN - MIN X TO ACCEPT IN PLOT
C       XMX - MAX X FOR PLOT
C       YMN - MIN Y TO ACCEPT IN PLOT
C       YMX - MAX Y FOR PLOT
C
C             IF XMN AND XMX OR YMN AND YMX ARE ZERO THEN THE
C             PLOT IS NOT FORCED, AND THE PACKAGE WILL FIND THE
C             RANGES FOR X AND/OR Y FROM THE STORED DATA.
C
C
      COMMON/PPBLK/ ITYPE,IWRITE,XMIN,XMAX,YMIN,YMAX,ICALL,IPOS,
     *              KJ,IFIN,IC(50),IND(50),INUM(50),FORCEX,FORCEY
      LOGICAL FORCEX,FORCEY
C
C
C   INITIALIZE
C
      IWRITE=IFLU
      ITYPE=NTYPE
      XMIN=XMN
      XMAX=XMX
      YMIN=YMN
      YMAX=YMX
      FORCEX=XMIN.NE.0. .OR. XMAX.NE.0.
      FORCEY=YMIN.NE.0. .OR. YMAX.NE.0.
      IF(.NOT. FORCEY) GO TO 1
      IF(ITYPE.EQ.1.OR.ITYPE.EQ.2.OR.ITYPE.EQ.5) GO TO 1
      IF(YMIN.GT.0. .AND. YMAX.GT.0.) GO TO 1
      WRITE(IWRITE,100)
      FORCEY=.FALSE.
1     IF(.NOT. FORCEX) GO TO 3
      IF(ITYPE.EQ.1.OR.ITYPE.EQ.3.OR.ITYPE.EQ.6) GO TO 3
      IF(XMIN.GT.0. .AND. XMAX.GT.0.) GO TO 3
      FORCEX=.FALSE.
      WRITE(IWRITE,100)
100   FORMAT(/' **WARNING** FORCE PARAMETERS .LE. 0, PLOT NOT FORCED')
3     ICALL=0
      IPOS=1
      KJ=1
      IFIN=-1
      RETURN
      END
      SUBROUTINE PPLOT(XLEN,PHITE,ITITX,ITITY,WORK)
C
C   ROUTINE THAT PERFORMS THE PLOTTING FOR THE PACKAGE
C
C
C      WHERE:
C
C      XLEN - THE PLOT LENGTH IN INCHES (MAX - 11)
C     PHITE - THE PLOT HEIGHT IN INCHES
C     ITITX - A VECTOR CONTAINING X AXIS ANNOTATION IN A1 FORMAT.
C             A ZERO (0) INDICATES END OF TITLE.
C     ITITY - A VECTOR CONTAINING Y AXIS ANNOTATION IN A1 FORMAT.
C             A ZERO (0) INDICATES END OF TITLE.
C      WORK - THE SAME WORK VECTOR USED IN THE PPSTOR STAGE.
C
C
      DIMENSION ITITX(1),ITITY(1),WORK(1),XLABEL(12)
      DIMENSION LINE(112),CODES(6)
      DIMENSION PLC(11),FORMAT(7)
      COMMON/PPBLK/ ITYPE,IWRITE,XMIN,XMAX,YMIN,YMAX,ICALL,IPOS,
     *              KJ,IFIN,IC(50),IND(50),INUM(50),FORCEX,FORCEY
      LOGICAL FORCEX,FORCEY
C
      DATA FORMAT/'(1X,','A1,F','XXXX',',1X,','1H+,','XXXX',
     + 'A1) '/
      DATA CODES/'6.5','6.4','6.3','6.2','6.1','6.0'/
      DATA PLC/'  12','  22','  32','  42','  52'
     * ,'  62','  72','  82','  92',' 102',' 112'/
      DATA MIX/'X'/,IMNS/'-'/,IOH/'O'/,IBLK/' '/,IPLS/'+'/
C  CLIP PLOT AT SINGLE PRECISION MAXIMUM
      DATA SPMAX/1.E38/
C
C  PLOT STORED VALUES
C
      IF(IPOS .GT. 1) GO TO 8
      WRITE(IWRITE,702)
702   FORMAT(/' ***ERROR*** NO PLOTS STORED')
      RETURN
C
C   GET X AND Y AXES SET UP
C
8     IF(.NOT. FORCEX) GO TO 81
      IF(ITYPE.NE.2.AND.ITYPE.NE.4.AND.ITYPE.NE.5.AND.ITYPE.NE.7)
     * GO TO 81
      XMIN=ALOG10(XMIN)
      XMAX=ALOG10(XMAX)
81    IF(.NOT. FORCEY) GO TO 80
      IF(ITYPE.NE.3.AND.ITYPE.NE.4.AND.ITYPE.NE.6.AND.ITYPE.NE.7)
     * GO TO 80
      YMIN=ALOG10(YMIN)
      YMAX=ALOG10(YMAX)
C   FIND X AND Y RANGE
80    IF(FORCEX .AND. FORCEY) GO TO 82
      IF(.NOT. FORCEY) YMAX=-SPMAX
      IF(.NOT. FORCEY) YMIN=SPMAX
      IF(.NOT. FORCEX) XMIN=SPMAX
      IF(.NOT. FORCEX) XMAX=-SPMAX
      DO 83 I=1,IFIN,2
      IF(.NOT. FORCEY) YMAX=AMAX1(YMAX,WORK(I+1))
      IF(.NOT. FORCEY) YMIN=AMIN1(YMIN,WORK(I+1))
      IF(.NOT. FORCEX) XMAX=AMAX1(XMAX,WORK(I))
      IF(.NOT. FORCEX) XMIN=AMIN1(XMIN,WORK(I))
83    CONTINUE
C
82    ILEN=XLEN+.49
      ILEN=MAX0(ILEN,1)
      ILEN=MIN0(ILEN,11)
      PLEN=ILEN*10
      XMX=XMAX
      XMN=XMIN
      YMX=YMAX
      YMN=YMIN
      IF(ITYPE.NE.3.AND.ITYPE.NE.4) GO TO 74
      IF(YMN .GE. 0.) YMIN=IFIX(YMN)
      IF(YMN .LT. 0.) YMIN=IFIX(YMN)-1
      IF(YMX .GE. 0.) YMAX=IFIX(YMX)+1
      IF(YMX .LT. 0.) YMAX=IFIX(YMX)
74    IF(ITYPE.NE.2.AND.ITYPE.NE.4) GO TO 84
      IF(XMN .GE. 0.) XMIN=IFIX(XMN)
      IF(XMN .LT. 0.) XMIN=IFIX(XMN)-1
      IF(XMX .GE. 0.) XMAX=IFIX(XMX)+1
      IF(XMX .LT. 0.) XMAX=IFIX(XMX)
84    A=YMAX-YMIN
      B=XMAX-XMIN
      IF(A .NE. 0.) GO TO 85
      YRANGE=.5
      YMIN=YMIN-1
      YMAX=YMAX+1
      GO TO 86
85    YRANGE=1./(YMAX-YMIN)
86    IF(B .NE. 0.) GO TO 87
      XRANGE=.5
      XMIN=XMIN-1
      XMAX=XMAX+1
      GO TO 88
87    XRANGE=1./(XMAX-XMIN)
C   COMPUTE INFO FOR AXES
88    XAXIS=(XMAX-XMIN)/FLOAT(ILEN)
      ILEN1=ILEN+1
      XLABEL(1)=XMIN
      DO 89 I=1,ILEN
89    XLABEL(I+1)=XLABEL(I)+XAXIS
      IF(ITYPE.NE.4.AND.ITYPE.NE.2.AND.ITYPE.NE.7.AND.ITYPE.NE.5)
     * GO TO 90
      DO 91 I=1,ILEN1
91    XLABEL(I)=10.**XLABEL(I)
90    XTOP=XRANGE*PLEN
      YLIM=6.*PHITE
      YLIM=IFIX(YLIM)
      YTOP=YLIM*YRANGE
      IYLIM=IFIX(YLIM)+1
      YAXIS=(YMAX-YMIN)/YLIM
      DO 92 I=1,IFIN,2
      WORK(I)=(WORK(I)-XMIN)*XTOP + 1.5
      WORK(I+1)=(WORK(I+1)-YMIN)*YTOP + 1.5
92    CONTINUE
C
C   START PLOTTING
C
      III=MIN0(6,MAX0(1,IFIX(ALOG10(YMAX)+2)))
      IF(ITYPE.EQ.3.OR.ITYPE.EQ.4.OR.ITYPE.EQ.6.OR.ITYPE.EQ.7)
     * III=MIN0(6,MAX0(1,IFIX(YMAX+2)))
      FORMAT(3)=CODES(III)
      FORMAT(6)=PLC(ILEN)
C   GET AXIS ANNOTATION INFO
      ILEN10=ILEN*10
      ICH=0
98    ICH=ICH+1
      IF(ITITX(ICH) .NE. 0) GO TO 98
      IXT=ICH-1
      ICH=0
97    ICH=ICH+1
      IF(ITITY(ICH) .NE. 0) GO TO 97
      IYT=ICH-1
      IF(IXT .EQ. 0) GO TO 95
      IXM=MIN0(IXT,ILEN10)
      ISP=(ILEN10-IXM+2)/2
      WRITE(IWRITE,803) (IBLK,I=1,ISP),(ITITX(I),I=1,IXM)
803   FORMAT(T9,122A1)
95    IYM=MIN0(IYT,IYLIM)
      ISY=(IYLIM-IYM)/2
C
C   PLOT X AXIS
C
      WRITE(IWRITE,801) (XLABEL(KKK),KKK=2,ILEN1,2)
801   FORMAT(T4,5(11X,1PE9.2))
      WRITE(IWRITE,800) (XLABEL(KKK),KKK=1,ILEN1,2)
800   FORMAT(T5,5(1PE9.2,11X),1PE9.2)
      DO 93 I=1,ILEN10
      LINE(I)=IMNS
      IF(MOD(I,10) .EQ. 0) LINE(I)=IPLS
93    CONTINUE
      ILEN11=ILEN10+1
      LINE(ILEN11)=IOH
      WRITE(IWRITE,802) (LINE(I),I=1,ILEN11)
802   FORMAT(T10,'O+',112A1)
C
C   DO ONE Y LINE AT A TIME
C
      ILEN12=ILEN11+1
      DO 100 III=1,IYLIM
      II=IYLIM-III+1
      DO 101 LL=1,ILEN11
101   LINE(LL)=IBLK
      M=II-1
      ITEST=MOD(M,6)
C   SCAN WORK ARRAY FOR VALUES OF Y EQUAL TO THIS Y VALUE
      DO 102 JJ=1,IFIN,2
      IY=WORK(JJ+1)
      IF(IY .NE. II) GO TO 102
      IX=WORK(JJ)
      JJ2=(JJ+1)/2
      DO 103 IP=2,IPOS
      IPS=IP-1
      IF(IND(IPS).LE.JJ2.AND.INUM(IPS).GE.JJ2) GO TO 104
103   CONTINUE
104   ICH=IC(IPS)
      IF(LINE(IX) .EQ. IBLK) GO TO 105
      IF(LINE(IX) .EQ. ICH) GO TO 102
      ICH=MIX
105   LINE(IX)=ICH
102   CONTINUE
      LINE(ILEN12)=IMNS
      ICH=IBLK
      IF(IYT.EQ.0 .OR. III.LE.ISY .OR. III.GT.IYT+ISY) GO TO 99
      ICH=ITITY(III-ISY)
99    IF(ITEST.EQ.0 .OR. III.EQ.1) GO TO 111
C   PRINT A LINE
      WRITE(IWRITE,900) ICH,(LINE(KKK),KKK=1,ILEN12)
900   FORMAT(1X,A1,T10,'-',112A1)
      GO TO 100
C   PRINT A LABELLED LINE
111   YY=YMIN+M*YAXIS
      IF(ITYPE.EQ.3.OR.ITYPE.EQ.4.OR.ITYPE.EQ.6.OR.ITYPE.EQ.7)
     * YY=10.**YY
      DO 110 IJ=1,ILEN11,10
      IF(LINE(IJ) .EQ. IBLK) LINE(IJ)=IPLS
110   CONTINUE
      LINE(ILEN12)=IPLS
      WRITE(IWRITE,FORMAT) ICH,YY,(LINE(KKK),KKK=1,ILEN12)
100   CONTINUE
C
C   PRINT X AXIS
C
      DO 96 I=1,ILEN10
      LINE(I)=IMNS
      IF(MOD(I,10) .EQ. 0) LINE(I)=IPLS
96    CONTINUE
      LINE(ILEN11)=IOH
      WRITE(IWRITE,802) (LINE(I),I=1,ILEN11)
      WRITE(IWRITE,800) (XLABEL(KKK),KKK=1,ILEN1,2)
      WRITE(IWRITE,801) (XLABEL(KKK),KKK=2,ILEN1,2)
      IF(IXT.NE.0) WRITE(6,803) (IBLK,I=1,ISP),(ITITX(I),I=1,IXM)
C
C  CLEAN UP PARAMETERS
C
      ICALL=0
      IPOS=1
      KJ=1
      IFIN=-1
      RETURN
      END
      SUBROUTINE PPSETU(IX,X,Y,ICAR,XLAB,YLAB,IPLOT)
      REAL X(100),Y(100),XLAB(80),YLAB(80)
      COMMON/UAB/WORK(3656)
      XH=0.
      XL=0.
      YL=0.
      YH=0.
      IF(IPLOT.EQ.-2)GO TO 200
      IF(IPLOT.EQ.0)GO TO 200
      IF(IPLOT.EQ.-1)GO TO 100
      YA=0.
      YL=Y(1)
      YH=Y(1)
      DO 10 I=1,IX
      IF(Y(I).LT.YL)YL=Y(I)
      IF(Y(I).GT.YH)YH=Y(I)
10    YA=YA+Y(I)
      YA=YA/IX
      IF(0.90*YA.LT.YL)YL=0.9*YA
      IF(1.10*YA.GT.YH)YH=1.10*YA
100   CALL PPINIT(1,6,XL,XH,YL,YH)
200   CALL PPSTOR(IX,X,Y,ICAR,WORK)
      IF(IPLOT.LT.0)GO TO 400
      WRITE(6,600)
600   FORMAT('1  ')
      CALL PPLOT(6.0,4.0,XLAB,YLAB,WORK)
400   CONTINUE
      RETURN
      END
      SUBROUTINE PPSTOR(NP,X,Y,ICHR,WORK)
C
C   STORAGE ROUTINE FOR PRINTER PLOTTING PACKAGE
C
C
C      WHERE:
C
C        NP - NUMBER OF POINTS TO STORE (FOR THIS PASS)
C         X - A VECTOR OF X VALUES
C         Y - A VECTOR OF Y VALUES
C      ICHR - THE PLOTTING CHARACTER
C      WORK - WORK VECTOR OF ATLEAST 2 TIMES THE TOTAL NUMBER
C                    OF POINTS STORED DURING THIS STAGE
C
C
C    PPSTOR MAY BE CALLED AS MANY AS 50 TIMES (TO OVERLAY
C    CURVES ON ONE ANOTHER) AFTER A CALL TO PPINIT.
C    A CALL TO PPLOT WILL PLOT THE DATA STORED BY THE PPSTOR
C    STAGE.
C
C
      DIMENSION X(1),Y(1),WORK(1)
      COMMON/PPBLK/ ITYPE,IWRITE,XMIN,XMAX,YMIN,YMAX,ICALL,IPOS,
     *              KJ,IFIN,IC(50),IND(50),INUM(50),FORCEX,FORCEY
      LOGICAL FORCEX,FORCEY
C
C   STORE THE VALUES FOR FUTURE PLOTTING
C
      N=NP
      IF(IPOS .EQ. 51) GO TO 60
      J=IFIN
      ISTART=ICALL
      DO 4 I=1,N
      IF(FORCEX .AND. (X(I).GT.XMAX .OR. X(I).LT.XMIN)) GO TO 4
      IF(FORCEY .AND. (Y(I).GT.YMAX .OR. Y(I).LT.YMIN)) GO TO 4
      J=J+2
      WORK(J)=X(I)
      WORK(J+1)=Y(I)
4     CONTINUE
      IF(J .EQ. IFIN) RETURN
      IFIN=J
      N=(J+1)/2-ICALL
C   CONVERT X AND Y VALUES IF NECESSARY
      GO TO(40,10,20,30,10,20,30),ITYPE
C
10    DO 11 IN=1,N
      I=2*(IN+ISTART)-1
      IF(WORK(I) .LE. 0.) GO TO 50
11    WORK(I)=ALOG10(WORK(I))
      GO TO 40
C
20    DO 21 IN=1,N
      I=2*(IN+ISTART)
      IF(WORK(I) .LE. 0.) GO TO 50
21    WORK(I)=ALOG10(WORK(I))
      GO TO 40
C
30    DO 31 IN=1,N
      I=2*(IN+ISTART)-1
      IF(WORK(I).LE.0. .OR. WORK(I+1).LE.0.) GO TO 50
      WORK(I)=ALOG10(WORK(I))
31    WORK(I+1)=ALOG10(WORK(I+1))
C
C   RECORD OTHER INFO
C
40    IC(IPOS)=ICHR
      IND(IPOS)=ISTART+1
      INUM(IPOS)=N+ISTART
      ICALL=ICALL+N
      IPOS=IPOS+1
      RETURN
C
C  ERROR MEAASGES
C
50    WRITE(IWRITE,700)
700   FORMAT(/' ***ERROR*** NEGATIVE OR ZERO VALUE ENTERED FOR '
     * ,'LOG PLOT.')
      GO TO 9
60    WRITE(IWRITE,701)
701   FORMAT(/' ***ERROR*** MAX PLOTS EXCEEDED - NO PLOT PRODUCED')
9     ICALL=0
      IPOS=1
      KJ=1
      IFIN=-1
      RETURN
      END
C     ****************************************************
C
C     THIS VERSION OF GRAPH3D USES UNITS 9, 10, AND 11
C     AS UNFORMATTED TEMPARARY FILES.
C
C     ***************************************************
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/PROP/GAMMA,W,G1
      COMMON/LSAVE/LPRI
      COMMON/SAVEU2/TU2
      COMMON/PASRAD/RADIUS(100,17,18)
      COMMON/PASS/K,L,R4
      COMMON/UAB/WORK(3696)
      REAL TU1(1728),TU2(1728)
      REAL JAC(50)
      REAL RTZ1(3,17,18), RTZ2(3,17,18),RTZG(20,20),X(20)
      REAL R(20),U1(5,17,18),U2(5,17,18),UT(5)
      REAL UA(120),RQ(20)
      REAL XPLOT(100),YPLOT(100)
      REAL XLAB(15),YLAB1(7),YLAB2(16)
      EQUIVALENCE (TU1,U1),(TU2,U2),(UA,RQ)
      DATA XLAB/'A','X','I','A','L',' ','D','I','S','T','A',
     1          'N','C','E',0/
      DATA YLAB1/'R','A','D','I','U','S',0/
      DATA YLAB2/'S','T','A','T','I','C',' ','P','R','E',
     1           'S','S','U','R','E',0/
      DATA ICARP,ICARS/'P','S'/
C
      CALL SETUP
      CALL OPEN
      IBUG=0
      LPRI=50
      NUMIT=3
      IR=1
      IW=2
C
C     WRITE GRID POINT OUTPUT - UNIT 7
C     L,K,J    =R,T,Z INDICES
C     R(L)     =R-COORDINATE
C     RTZG(K,L)=T-COORDINATE
C     X(L)     =Z-COORDINATE
C     U1(I,K,L)=RHO*R*(1,UR,UTH,UZ),R*ET
C     U2(1,K,L)=PT,MR,MTH,MZ,TT
C
      CALL SETCOR(IBUG,X,R,RTZG)
      DO 4000 J=1,NX
      CALL IO(IR,UA,0,J,TU1,JAC)
      READ(11)X,R,RTZG
      DO 3500 L=1,NR
      DO 3500 K=1,NTH
      RADIUS(J,K,L)=R(L)
3500  CONTINUE
      CALL PTASS(0,-5,J,J,1,NR,0,0)
      DO 4000 L=1,NR
      DO 4000 K=1,NTH
C
      WRITE(7)L,K,J,R(L),RTZG(K,L),X(L),(U1(I,K,L),I=1,5),
     1(U2(I,K,L),I=1,5)
C
4000  CONTINUE
      REWIND 11
      WRITE(6,6020)
6020  FORMAT(1X,' GRAPH3D OUTPUT AT GRID POINTS')
      NC=NR/2+1
      CALL PTASS(0,5,N1,N1,NC,NC,0,0)
      CALL PTASS(0,5,N2,N2,NC,NC,0,0)
C
C     CALCULATE STREAMLINE POSITIONS
C
      NR1=NR-1
      NTH1=NTH-1
      NA=N1+1
      CALL SETRTZ(RTZ1,X,R,RTZG)
      CALL IO(IR,UA,0,N1,TU1,JAC)
      CALL PPROP(U1,RQ)
      KEY=1
      CALL CPROP(U1,RTZ1,RTZ2,KEY)
      WRITE(10)RTZ2
      CALL IO(IW,UA,0,N1,TU1,JAC)
      DO 5 L=1,NR
      DO 5 K=1,NTH
5     RADIUS(N1,K,L)=R(L)
C
      DO 900 J=NA,N2
      CALL PPROP(U1,RQ)
      CALL IO(IR,UA,0,J,TU2,JAC)
      CALL SETRTZ(RTZ2,X,R,RTZG)
      CALL PPROP(U2,RQ)
      DO 6 L=1,NR
      DO 6 K=1,NTH
6     RADIUS(J,K,L)=RTZ2(1,K,L)
C
      DO 100 L=1,NR
      DO 100 K=1,NTH
      R1=RTZ1(1,K,L)
      T1=RTZ1(2,K,L)
      X1=RTZ1(3,K,L)
      X2=RTZ2(3,K,L)
      DELX=X2-X1
      UR=U1(2,K,L)/U1(4,K,L)
      UTH=U1(3,K,L)/U1(4,K,L)
      R2=R1+DELX*UR
      T2=T1+DELX*UTH
      R2=RTZ2(1,K,L)
      T2=RTZ2(2,K,L)
      DO 10 IP=1,NUMIT
      KP=K
      LP=L
      CALL SMLINT(U2,RTZG,R,X,R2,T2,UT,X2,IBUG,KP,LP)
      DELX=X2-X1
      R2=R1+DELX/2.*(UR+UT(2)/UT(4))
      T2=T1+DELX/2.*(UTH/R1+UT(3)/UT(4)/R2)
10    CONTINUE
      DO 15 I=1,5
15    U1(I,K,L)=UT(I)
      RADIUS(J,K,L)=R2
      RTZ1(1,K,L)=R2
      RTZ1(2,K,L)=T2
      RTZ1(3,K,L)=X2
100   CONTINUE
      KEY=1
      CALL CPROP(U1,RTZ1,RTZ2,KEY)
      WRITE(10)RTZ2
      CALL IO(IW,UA,0,J,TU1,JAC)
      L=LPRI
900   CONTINUE
C
      WRITE(6,6030)
6030  FORMAT(1H1,' GRAPH3D OUTPUT ON STREAM SURFACES')
      CALL PTASS(0,5,N1,N1,NC,NC,0,0)
      CALL PTASS(0,5,N2,N2,NC,NC,0,0)
C
C     WRITE STREAMLINE OUTPUT - UNIT 8
C     STREAMLINES START AT L.E. GRID POINTS
C     L,K,J,U1,U2=AS ABOVE
C     RTZ1(1,K,L)=R-COORDINATE
C     RTZ1(2,K,L)=T-COORDINATE
C     RTZ1(3,K,L)=Z-COORDINATE
C
      REWIND 10
      RH=0.
      RL=1000.
      XH=-1000.
      XL=1000.
      DO 920 J=N1,N2
      CALL IO(IR,UA,0,J,TU1,JAC)
      READ(10)RTZ1
      CALL PTASS(0,-5,J,J,1,NR,0,0)
      DO 4500 L=1,NR
      DO 4500 K=1,NTH
C
      WRITE(8)L,K,J,RTZ1(1,K,L),RTZ1(2,K,L),RTZ1(3,K,L),
     1(U1(I,K,L),I=1,5),(U2(I,K,L),I=1,5)
C
4500  CONTINUE
      DO 910 K=1,NTH
      DO 910 L=1,NR
      X1=RTZ1(3,K,L)
      R1=RTZ1(1,K,L)
      R4=R1
      CALL FINDP(U1,PST,VV)
      U1(1,K,L)=PST
      U1(2,K,L)=X1
      IF(X1.LT.XL)XL=X1
      IF(X1.GT.XH)XH=X1
      IF(R1.LT.RL)RL=R1
      IF(R1.GT.RH)RH=R1
910   CONTINUE
      CALL IO(IW,UA,0,J,TU1,JAC)
920   CONTINUE
C
C     MAXE STREAMLINE PLOTS ON S2 SURFACES
C
      CALL PPINIT(1,6,XL,XH,RL,RH)
      DO 950 K=1,NTH,NTH1
      WRITE(6,6000)K
6000  FORMAT('1             S2 SURFACE NUMBER',I5)
      REWIND 10
      DO 948 J=N1,N2
      J1=0
      READ(10)RTZ1
      DO 945 LQQ=1,NR
      L=LQQ
      IF(LQQ.GE.NR)L=NR
      J1=J1+1
      XPLOT(J1)=RTZ1(3,K,L)
      YPLOT(J1)=RTZ1(1,K,L)
945   CONTINUE
      CALL PPSTOR(J1,XPLOT,YPLOT,'*',WORK)
948   CONTINUE
      CALL PPLOT(6.0,7.0,XLAB,YLAB1,WORK)
950   CONTINUE
C
C     MAKE STATIC PRESSURE PLOTS ON S1 SURFACES
C
      DO 990 L=2,NR1
      WRITE(6,6002)L
6002  FORMAT('1 '//'      STATIC PRESSURE DIVIDED BY PREF',
     1            /'      FOR S1 SURFACE ',I5)
      DO 970 K=1,NTH,NTH1
      PMAX=0.
      PMIN=0.
      DO 960 J=N1,N2
      J1=J-N1+1
      CALL IO(IR,UA,0,J,TU1,JAC)
      PST=U1(1,K,L)*GAMMA
      IF(PST.GT.PMAX)PMAX=PST
      XPLOT(J1)=U1(2,K,L)
      YPLOT(J1)=PST
960   CONTINUE
      ICAR=ICARP
      IF(K.GT.1)ICAR=ICARS
      IF(K.GT.1)GO TO 965
      CALL PPINIT(1,6,XL,XH,PMIN,PMAX)
965   CALL PPSTOR(J1,XPLOT,YPLOT,ICAR,WORK)
970   CONTINUE
      CALL PPLOT(6.0,7.0,XLAB,YLAB2,WORK)
990   CONTINUE
C
      STOP
      END
      SUBROUTINE CPROP(U,R,R3,KEY)
      REAL U(5,17,18),R(3,17,18)
      REAL R3(3,17,18)
      COMMON/GRID/NX,NTH,NR
      COMMON/PROP/GAMMA,W,G1
      DO 50 L=1,NR
      DO 50 K=1,NTH
      R1=R(1,K,L)
      RHOR=U(1,K,L)*R1
      IF(KEY.EQ.1)U(3,K,L)=U(3,K,L)+W*R1
      DO 10 I=2,5
10    U(I,K,L)=U(I,K,L)*RHOR
      U(1,K,L)=RHOR
30    CONTINUE
      R3(1,K,L)=R1
      R3(3,K,L)=R(3,K,L)
      R3(2,K,L)=R(2,K,L)*R1
50    CONTINUE
      RETURN
      END
      SUBROUTINE FINDP(U,P,V)
C       THIS ROUTINE CALCULATES PRESSURE FROM CONSERVATION VARIABLES
      REAL U(5,17,18)
      COMMON /PROP/GAMMA,WSAV,G1
      COMMON/ROTOR/TITLE(80),NBLADE,ICONS,IUSE,IROT
      COMMON /PASS/K,L,R
C       CALCULATE  SQUARED VELOCITY V=(UR**Z + UTH**2 + UZ**2)*RHO*R
      V = (U(2,K,L)**2+U(3,K,L)**2+U(4,K,L)**2)/U(1,K,L)
      IF(IROT.NE.0) GO TO 10
      P = (U(5,K,L)-0.5*V)/R*G1
      RETURN
C
C       CALCULATE P USING CONSTANT RHOTHAPLY
C       FAR UPSTREAM ENTHALPY IS 1/(GAMMA-1)
C
10    CONTINUE
      RHOR=U(1,K,L)
      RHIN = 1.0/RHOR
      UTH=U(3,K,L)*RHIN
      U2=V*0.5*RHIN
      E=(1./G1+R*WSAV*UTH-U2)/GAMMA
      P=RHOR/R*E*G1
      U(5,K,L)=RHOR*(E+U2)
      RETURN
      END
      SUBROUTINE IO(IRW,U,IUORUT,J,V,JAC)
C
C     THIS ROUTINE DOES ALL DISK I/O FOR PROBLEM MATRIX
C     IRW=1 => READ   IRW= 2 => WRITE
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /US1/U1(1728,100),UA1(120,100),JAC1(20,100)
      REAL JAC1
      REAL V(1728),JAC(50)
      REAL U(120)
      IF(IRW.EQ.2) GO TO 100
      DO  10 I=1,1728
10    V(I)=U1(I,J)
      DO  20 I=1,120
20    U(I)=UA1(I,J)
      DO  30 I=1,NR
30    JAC(I)=JAC1(I,J)
C
      RETURN
100   CONTINUE
      DO  200 I=1,1728
200   U1(I,J)=V(I)
      RETURN
      END
      SUBROUTINE OPEN
      COMMON /GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON /US1/U1(1728,100),UA1(120,100),JAC1(20,100)
      COMMON/TGSAV/TTSAV(18,18,100),TRSAV(18,18,100),TXSAV(18,18,100),
     *   CURSAV(18,2,2,100),COSNSV(18,3,2,100)
      REAL TU(1728),UA(150),JAC(50)
      REAL JAC1
      INTEGER XXS,XXE,XRS,XRE,RXS,RXE,RRS,RRE
      DATA XXS,XRS,RXS,RRS/19,37,55,73/
      XXE=NR+18
      XRE=NR+36
      RXE=NR+54
      RRE=NR+72
      DO 100 J=1,NX
      READ(2)((TTSAV(K,L,J),K=1,NTH),L=1,NR),
     *((TRSAV(K,L,J),TXSAV(K,L,J),K=1,NTH),L=1,NR),
     *(((CURSAV(L,KX,KY,J),KX=1,2),(COSNSV(L,KX,KY,J),KX=1,3),
     *KY=1,2),L=1,NR)
      READ(3)(UA1(I,J),I=1,NR),(UA1(I,J),I=XXS,XXE),
     *(UA1(I,J),I=XRS,XRE),(UA1(I,J),I=RXS,RXE),
     *(UA1(I,J),I=RRS,RRE)
      READ(3)(UA1(I,J),I=91,105)
      READ(3)(JAC1(L,J),L=1,NR)
100   CONTINUE
C.... READ OLD SOLUTION.
      READ(1) DUM,DUM,DUM,DUM,DUM,DUM
      DO 10 J=1,NX
      READ(1)TU
      DO 5 I=1,1728
5     U1(I,J)=TU(I)
10    CONTINUE
      READ(1,END=11)IFLIPX,IFLIPT,IFLIPR,ICYCLE
11    CONTINUE
      REWIND 1
      REWIND 2
      REWIND 3
      RETURN
      END
      SUBROUTINE PPROP(U,R)
      REAL U(5,17,18),R(20)
      COMMON/GRID/NX,NTH,NR
      COMMON/PROP/GAMMA,W,G1
      COMMON/PASS/KK,LL,RR
      NTH1=NTH-1
      NR1=NR-1
C
C               CORRECT DUMMY POINTS
C
      DO 9 L=1,NR
      DO 9 K=1,NTH
      KK=K
      LL=L
      IF(L.EQ.1)LL=2
      IF(L.EQ.NR)LL=NR-1
      IF(K.EQ.1)KK=2
      IF(K.EQ.NTH)KK=NTH-1
      IF(L.EQ.1 .OR. L.EQ.NR)GO TO 2
      IF(K.GT.1 .AND. K.LT.NTH) GO TO 9
2     CONTINUE
      RR=R(L)
      CALL FINDP(U,P,V)
      RHOHO=U(5,KK,LL)/RR+P
      RHOIO=RHOHO-W*U(3,KK,LL)
      RHO=U(1,KK,LL)/RR
      ROTH=RHOIO/RHO
      LL=L
      KK=K
      RR=R(L)
      CALL FINDP(U,P,V)
      HO=ROTH+W*RR*U(3,K,L)/U(1,K,L)
      H=HO-0.5*V/U(1,K,L)
      E=H/GAMMA
      RHO=P/G1/E
      RHOR=RHO*RR
      ET=(E+0.5*V/U(1,K,L))*RHOR
      U(5,K,L)=ET
      DO 8 I=2,4
8     U(I,K,L)=U(I,K,L)/U(1,K,L)*RHOR
      U(1,K,L)=RHOR
9     CONTINUE
      DO 50 L=1,NR
      R1=R(L)
      DO 20 K=1,NTH
      RHOR=U(1,K,L)
      DO 10 I=2,5
10    U(I,K,L)=U(I,K,L)/RHOR
      U(3,K,L)=U(3,K,L)-W*R1
      UX=U(4,K,L)
C      UX=ABS(UX)
C      IF(UX.LT.0.1)UX=0.1
      U(4,K,L)=UX
C     IF(L.EQ.1)R1=R(2)
C     IF(L.EQ.NR)R1=R(NR1)
      U(1,K,L)=U(1,K,L)/R1
20    CONTINUE
50    CONTINUE
      DO 60 K=1,NTH,NTH1
      K1=2
      IF(K.NE.1)K1=NTH1
      DO 55 I=1,5
      IF(K.EQ.1)U(I,K,1)=U(I,2,1)-(U(I,3,1)-U(I,2,1))
      IF(K.EQ.1)U(I,K,NR)=U(I,2,NR)-(U(I,3,NR)-U(I,2,NR))
      IF(K.EQ.NTH)U(I,K,1)=U(I,K,1)-(U(I,K-1,1)-U(I,K-2,1))
      IF(K.EQ.NTH)U(I,K,NR)=U(I,K,NR)-(U(I,K-1,NR)-U(I,K-2,NR))
55    CONTINUE
      DO 56 L=1,NR
      DO 56 I=1,5
56    U(I,K,L)=(U(I,K,L)+U(I,K1,L))*0.5
60    CONTINUE
      DO 70 K=1,NTH
      DO 70 I=1,5
      U(I,K,1)=0.5*(U(I,K,1)+U(I,K,2))
      U(I,K,NR)=0.5*(U(I,K,NR)+U(I,K,NR-1))
70    CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(U,IMOD,MSG,J,LS,LSS,DMD,R1)
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/PROP/GAMMA,W,G1
      COMMON/PASRAD/RADIUS(100,17,18)
      REAL U(5,17,18),MSG(5,5)
      REAL R1(18)
      REAL R(18)
      REAL AV(5)
      TODEG=57.29578
      G2=G1/GAMMA
      NR1=NR-1
      NTH1=NTH-1
      NR2=NR-2
      RTIP=(R1(NR)+R1(NR-1))*0.5
      RHUB=(R1(1)+R1(2))*0.5
      IF(R1(1).EQ.R1(2))RHUB=0.
      DELR=RTIP-RHUB
      DO 40 L=LS,LSS
      WRITE(6,300)J,L,(K,K=1,NTH)
      SP=(R1(L)-RHUB)/DELR
      DO 10 I=1,5
      WRITE(6,301)MSG(I,IMOD),(U(I,K,L),K=1,NTH)
      SP=0.
      AV1=0.
      DO 5 K=1,NTH
      AV1=AV1+U(I,K,L)
5     SP=SP+RADIUS(J,K,L)
      AV(I)=AV1/NTH
      SP=SP/NTH
10    SP=(SP-RHUB)/DELR
      WRITE(6,303)(MSG(I,IMOD),I=1,5),SP,(AV(I),I=1,5)
40    CONTINUE
300   FORMAT(///I4,I3,3X,17(5X,I2))
301   FORMAT(5X,A4,1X,17F7.3)
303   FORMAT(/,T2,'THETA AVERAGE PROPERTIES.',/,
     1T6,'%SPAN',5A7/
     2T5,6(F6.3,1X))
      WRITE(6,305)
305   FORMAT(/,T3,'L  RADIUS MACH NUMBERS         P  ON  T  ON',
     1'  ANGLES',/,T13,
     2'MERID  TANGEN TOTAL  P REF  T REF  WHIRL  MERID')
306   FORMAT(T2,I2,2(2X,F5.3),2X,F5.2,3(2X,F5.3),2(2X,F5.1),2X,F5.3)
C
      DO 90 L=1,NR
      R(L)=R1(L)
      DO 50 I=1,5
      AV1=0.
      DO 45 K=1,NTH
45    AV1=AV1+U(I,K,L)
      AV(I)=AV1/NTH
50    CONTINUE
      TMR=(AV(2)**2 + AV(4)**2)
      TMR=SQRT(TMR)
      TM=TMR**2 +AV(3)**2
      TM=SQRT(TM)
      AV(1)=AV(1)*GAMMA
      ATW=ATAN(AV(3)/TMR)*TODEG
      ATR=ATAN(AV(2)/AV(4))*TODEG
      WRITE(6,306)L,R(L),TMR,AV(3),TM,AV(1),AV(5),ATW,ATR
90    CONTINUE
C
      RETURN
      END
      SUBROUTINE PTASS(IPRT,NM,JSS1,JSS2,LQ1,LQ2,ISEC,KSAV)
C
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/PROP/GAMMA,W,G1
      COMMON/ROTOR/TITLE(80),NBLADE,ICONS,IUSE,IROT
      COMMON/PASS/K,L,R
      COMMON/PROUT/PRPS(100),PRSS(100)
      COMMON/PASRAD/RADIUS(100,17,18)
      COMMON/SAVEU2/TU2(1728)
      REAL UA(120),U(5,17,18),R1(18),MDOT
      REAL XR1(18),XX1(18),RX1(18),RR1(18)
      REAL JAC(50)
      REAL MSG(5,5),MX,M,SYS(5)
      EQUIVALENCE (U,TU2)
      EQUIVALENCE (UA(1) , R1(1)),(UA(19),XX1(1))
      EQUIVALENCE (UA(37),XR1(1)),(UA(55),RX1(1))
      EQUIVALENCE (UA(73),RR1(1))
      DATA MSG/4HRHOP,4HURP ,4HUTHP,4HUXP ,4HETP ,
     *         4HRHO ,4HUR  ,4HUTH ,4HUX  ,4HET  ,
     *         4HRHO ,4HMR  ,4HMTH ,4HMX  ,4HP   ,
     *         4HPT  ,4HMR  ,4HMTH ,4HMX  ,4HTT  ,
     *         4HPT  ,4HMR  ,4HMTH ,4HMX  ,4HTT  /
      DATA SYS/4HCONS,4HABS ,4HREL ,4HREL ,4HABS /
C
C     MODE=1  =CONSERVATION VARIABLES
C     MODE=2  = PHYSICAL VARIABLES IN ABS SYSTEM
C     MODE=3  = RHO,MR,MTH,MX,P    IN REL SYSTEM
C     MODE=4  =PT,MR,MTK,MX,TT     IN REL SYSTEM
C     MODE=5  =PT,MR,MTH,MX,TT     IN ABS SYSTEM
C
      CALL IO(1,UA,0,1,TU2,JAC)
      KCY=TU2(1728)
      G2=0.5*G1
      G3=GAMMA/G1
      MODE = IABS(NM)
C
      IF(NM.LT.0) GO TO 55
      WRITE(6,301)JSS1,JSS2,SYS(MODE)
301   FORMAT(///1X,T2,'GRID OUTPUT FOR',I4,' TO',I4,
     *' IN ',A4,' SYSTEM')
      WRITE(6,304)(TITLE(I),I=1,80)
304   FORMAT(1X,T1,80A1)
55    IMOD = MODE
      DO 500 J=JSS1,JSS2
      JJ=J
      IA=0
      CALL IO(1,UA,IA,JJ,TU2,JAC)
      JKK=3
      IF(J.EQ.1 .OR. J.EQ.NX)JKK=2
      L=NR/2
      K=2
      R=RADIUS(J,K,L)
      CALL FINDP(U,P,V)
      PRPS(J)=P
      K=NTH-1
      R=RADIUS(J,K,L)
      CALL FINDP(U,P,V)
      PRSS(J)=P
      DO 90 L=1,NR
      R=R1(L)
      DO 90 K=1,NTH
      CALL FINDP(U,P,V)
      E=(U(5,K,L)-0.5*V)/U(1,K,L)
      T=GAMMA*G1*E
      IF(T.GT.0.)GO TO 936
      WRITE(6,510)J,K,L
510   FORMAT(' PRESSURE LESS THAN ZERO AT ',3I5)
      T=1.0
936   CONTINUE
      C=SQRT(T)
      RHOR=U(1,K,L)
      U(4,K,L)=U(4,K,L)/C/RHOR
      U(2,K,L)=U(2,K,L)/C/RHOR
      U(3,K,L)=U(3,K,L)/C/RHOR
      U(5,K,L)=P
      U(1,K,L)=RHOR/R
90    CONTINUE
      DO 100 L=1,NR
      DO 100 K=1,NTH
      M=0.
      DO 91 I=2,4
      M=M+U(I,K,L)**2
91    CONTINUE
      T1=1.0+G2*M
      PT=U(5,K,L)*T1**G3
      E=U(5,K,L)/U(1,K,L)/G1
      T=E*GAMMA*G1
      TT=T*T1
      U(5,K,L)=TT
      U(1,K,L)=PT
100   CONTINUE
400   CONTINUE
      IF(NM.LT.0)GO TO 500
      CALL PRINT(U,IMOD,MSG,J,LQ1,LQ2,MDOT,R1)
500   CONTINUE
      RETURN
      END
      SUBROUTINE SETRTZ(RTZ,X,R,RTZG)
      COMMON/GRID/NX,NTH,NR
      REAL RTZ(3,17,18),X(20),R(20),RTZG(20,20)
      READ(9)X,R,RTZG
      NR1=NR-1
      DO 100 L=1,NR
      R1=R(L)
      X1=X(L)
      DO 40 K=1,NTH
      RTZ(1,K,L)=R1
      RTZ(2,K,L)=RTZG(K,L)
40    RTZ(3,K,L)=X1
100   CONTINUE
      RETURN
      END
      SUBROUTINE SETCOR(IBUG,X,R,RTZG)
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/LSAVE/LPRI
      REAL X(20),R(20),TP(20),TS(20),RTZG(20,20)
      REAL XSAVE(100,20),RSAVE(100,20),TPSAVE(100,20),TSSAVE(100,20)
      COMMON/SVPASS/XSAVE,RSAVE,TPSAVE,TSSAVE
      COMMON/GRID/NX,NTH,NR,DUM(7)
      COMMON/ROTOR/TITLE(80),NBLADE,ICONS,IUSE,IROT
      NTH1=NTH-1
      IB=0
      IF(IBUG.EQ.-1.OR.IBUG.EQ.1)IB=1
      DO 100 J=1,NX
      DO 10 L=1,NR
      X(L)=XSAVE(J,L)
      R(L)=RSAVE(J,L)
      TS(L)=TSSAVE(J,L)
      TP(L)=TPSAVE(J,L)
10    CONTINUE
      TS(1)=0.5*(TS(1)+TS(2))
      TP(1)=0.5*(TP(1)+TP(2))
      TS(NR)=0.5*(TS(NR)+TS(NR-1))
      TP(NR)=0.5*(TP(NR)+TP(NR-1))
      DO 50 L=1,NR
      RP=R(L)
      RP=1.
      TP(L)=TP(L)*RP
      TS(L)=TS(L)*RP
      IF(IB.EQ.1.AND.L.EQ.LPRI)WRITE(6,600)J,L,TP(L),TS(L)
      RP=1
      DELT=(TS(L)-TP(L))*RP/(NTH-2)
      DELT2=DELT/2.
      RTZG(1,L)=TP(L)*RP
      RTZG(NTH,L)=TS(L)*RP
      DO 20 K=2,NTH1
      RTZG(K,L)=RP*TP(L)+DELT2+DELT*(K-2)
      IF(IB.EQ.1.AND.L.EQ.LPRI)WRITE(6,600)L,K,RTZG(K,L)
600   FORMAT(2I5,2F10.5)
20    CONTINUE
50    CONTINUE
      R(1)=0.5*(R(1)+R(2))
      R(NR)=0.5*(R(NR)+R(NR-1))
      X(1)=0.5*(X(1)+X(2))
      X(NR)=0.5*(X(NR)+X(NR-1))
      WRITE(11)X,R,RTZG
      IF(J.LT.N1)GO TO 100
      IF(J.GT.N2)GO TO 100
      WRITE(9)X,R,RTZG
100   CONTINUE
      WRITE(9)NBLADE
      REWIND 9
      REWIND 11
      RETURN
      END
      SUBROUTINE SETUP
      REAL TITL(20),X(100),Y(100),X1(20),R1(20)
      REAL XSAVE(100,20),RSAVE(100,20),TPSAVE(100,20),TSSAVE(100,20)
      COMMON/SVPASS/XSAVE,RSAVE,TPSAVE,TSSAVE
      COMMON/GRID/NX,NTH,NR,DX,DTH,DR,DT,TAUX,TAUT,TAUR
      COMMON/BLADE/N1,N2,KUTTA,IEXPLE,IEXPTE
      COMMON/UPBOND/AZERO,RHZERO,PDOWN,UPMACH
      COMMON/PROP/GAMMA,W,G1
      COMMON/ROTOR/TITLE(80),NBLADE,IFCL,IUSE,IROT
      COMMON/PRESS/P(3,18,4),LR,IDLOW,IDHIGH
C
      IUSE=0
      TAUX=0.0
      TAUT=0.0
      TAUR=0.0
      READ(5,555)(TITLE(I),I=1,80)
555   FORMAT(80A1)
      READ(5,501)NBLADE
501   FORMAT(5I5)
      IF(NBLADE.LE.0)NBLADE=1
      READ(5,501)NX,NTH,NR
      READ(5,501)N1,N2,KUTTA,IEXPLE,IEXPTE
      READ(5,501)IDLOW,IDHIGH,IFCL,IROT
      IF(IDLOW.LE.0)IDLOW=1
      IF(IDHIGH.LE.0)IDHIGH=NR
      DX=2.0/(N2-N1)
      DR=1.0/(NR-2)
      DTH=1.0/(NTH-2)
      READ(5,502)GAMMA,W,CAPPA,PDOWN
502   FORMAT(8F10.5)
C
      G1=GAMMA-1
      NTH1=NTH-1
      NR1=NR-1
C
      READ(4) NDM,NXS,NTHS,NRS,NDM,NDM,NDM,IEXS,NDM,NRHS,
     * NRTS,NZLES,NZTES,NZLS,NZUS
      READ(4) DUM
      READ(4) (DUM,I=1,NRHS)
      READ(4) (DUM,I=1,NRHS)
      READ(4) (DUM,I=1,NRTS)
      READ(4) (DUM,I=1,NRTS)
      READ(4) (DUM,I=1,NZLES)
      READ(4) (DUM,I=1,NZLES)
      READ(4) (DUM,I=1,NZTES)
      READ(4) (DUM,I=1,NZTES)
      IF(IEXS.LE.0) GO TO 101
      READ(4) (DUM,I=1,NZLS)
      READ(4) (DUM,I=1,NZLS)
      READ(4) (DUM,I=1,NZUS)
      READ(4) (DUM,I=1,NZUS)
101   READ(4) ((XSAVE(J,L),J=1,NXS),L=1,NRS)
      READ(4) ((RSAVE(J,L),J=1,NXS),L=1,NRS)
      READ(4) ((TSSAVE(J,L),J=1,NXS),L=1,NRS)
      READ(4) ((TPSAVE(J,L),J=1,NXS),L=1,NRS)
      REWIND 4
      RETURN
      END
      SUBROUTINE SMLINT(U,RTZG,R,X,R1,T1,U2,X2,IBUG,KP,LP)
      COMMON/LSAVE/LPRI
      REAL U(5,17,18),RTZG(20,20),R(20),X(20),U2(5)
      REAL T(20)
      COMMON /GRID/NX,NTH,NR
      IB=0
      IF(IBUG.EQ.-1.OR.IBUG.EQ.2)IB=1
      NTH1=NTH-1
      NTH2=NTH-2
      NR1=NR-1
      LS=LP
      IF(LP.EQ.NR)GO TO 10
      DO 5 L=1,NR1
      RS=R(L)
      IF(R1.GT.RS)GO TO 5
      LS=L
      GO TO 10
5     CONTINUE
      LS=NR
10    CONTINUE
      LS1=LS-1
      IF(LP.EQ.1)LS1=1
      IF(LP.EQ.1)LS=2
      RP=R(LS)-R(LS1)
      IF(RP.EQ.0.)WRITE(6,5000)LS,LS1,KP,LP
5000  FORMAT('  LS LS1 KP LP',4I5)
      RP=(R1-R(LS1))/RP
      IF(LP.EQ.1)RP=0.
      IF(LP.EQ.NR)RP=1.
      TA=(RTZG(1,LS)-RTZG(1,LS1))*RP+RTZG(1,LS1)
      TB=(RTZG(NTH,LS)-RTZG(NTH,LS1))*RP+RTZG(NTH,LS1)
      DELTH=(TB-TA)/NTH2
      DEL2=DELTH/2
      DO 20 K=2,NTH1
20    T(K)=TA+DELTH*(K-2)+DEL2
C     T(1)=TA-DEL2
      T(1)=TA
C     T(NTH)=TB+DEL2
      T(NTH)=TB
      IF(KP.EQ.1)GO TO 205
      IF(KP.LT.NTH)GO TO 203
      IK=NTH
      TP=T(NTH)-T(NTH1)
      GO TO 30
203   CONTINUE
      IF(T1.GT.TA)GO TO 210
205   CONTINUE
      IK=2
      TP=0.
      GO TO 30
210   CONTINUE
      DO 25 K=2,NTH
      TQ=T(K)
      IF(T1.GT.TQ)GO TO 25
      IK=K
      TP=T1-T(IK-1)
      GO TO 30
25    CONTINUE
      IK=NTH
      TP=TB-T(NTH1)
30    CONTINUE
      IK1=IK-1
      TP=TP/(T(IK)-T(IK1))
      IF(IB.EQ.1)WRITE(6,507)T1,TP,TA,TB,KP,LP
507   FORMAT(' T1,TP,T',4F10.5,2I5)
      X2=(X(LS)-X(LS1))*RP+X(LS1)
      IF(IB.EQ.1)WRITE(6,508)RP,LS,TP,IK
508   FORMAT(' RP,LS,TP,IK',F10.5,I5,F10.5,I5)
      DO 50 I=1,5
      A1=U(I,IK1,LS1)
      A1=(U(I,IK1,LS)-A1)*RP+A1
      B1=U(I,IK,LS1)
      B1=(U(I,IK,LS)-B1)*RP+B1
      C1=(B1-A1)*TP+A1
      U2(I)=C1
50    CONTINUE
      RETURN
      END
